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Abstract 

In this work we present a calculation of the hamiltonian variables solv- 
ing the molecular dynamics equations of motion for a system of nuclear 
matter kept at fixed temperature by using the Nose-Hoover Thermostat 
and interacting via a semiclassical potential depending on both positions 
and momenta. 



1 Introduction 

At densities just below nuclear saturation density, there may be spatial struc- 
tures different from the uniformly distributed matter. In events like a supernova 
core collapse large quantities of neutrinos are produced. As they stream out of 
the star, coherent scattering out of spatial density fluctuations of neutron rich 
matter, known as pasta phases, can happen. These structures arise due to 
the competition among intermediate-range attractive and long-range repulsive 
forces p. Possible shapes include round nuclei, flat plates, rods, and spherical 
voids |2] ■ 

In this work we are interested in solving the dynamical equations for a nu- 
clear system interacting via a realistic potential model depending not only on 
positions, as usual, but on both positions and momenta and keeping the temper- 
ature fixed using the Nose-Hoover thermostatic]. We mimic the Pauli principle 
in a fermionic nucler system by adding a momentum dependent term to the 
potential py. In this way the characteristic phase space repulsion for fermionic 
nucleons (protons and neutrons) can be included, restricting some of the avail- 
able dynamical states for the individual particles. Solving this system allows 
for realistic configurations of neutron rich plasmas in the NVT ensemble. 
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2 Model Hamiltonian and Method 



We model a charge-neutral system with a fixed number of nucleons, A, and 
electrons. The electrons provide a neutralizing background and are described 
as an almost degenerate free Fermi gas, see below. The Hamiltonian under the 
Nose-Hoover method for the extended system in this case can be written as |3] 
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where V(Rij,Pij) is the potential which depends on both positions and mo- 
menta, described below, s is the extended position variable, p s is the momentum 
conjugate to s, Q is the thermal incrtial parameter corresponding to a coupling 
constant between the system and the thermostat taking a value Q ~ 10 6 -10 8 
MeV (/m/c) 2 , we use g = 3 A as a condition for generating the canonical en- 
semble in the classical molecular dynamics simulations, £ is the thermodynamic 
friction coefficient and j3 is defined as (3 = 1/k-oT. 

The total potential energy of the system, V, consists of a sum of two-body 
interactions 

V = Vhad + Vcoulomb + Vp a uli (2) 
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Here the distance between the particles in phase space is denoted by Rij = 
|Rj — Rj|, Pij — |Pj— Pj| and represents the ith-nucleon isospin projection 
on z-axis (r =+1 for protons and r = — 1 for neutrons). Vcoulomb corresponds 
to the screened Coulomb interaction. The screening length, A, that results 
from the slight polarization of the electron gas is arbitrarily set to A = 10 fm 
as in previous works 0. Vp au u is the Pauli potential that incorporates 

phase space repulsion for fermions by means of the Kronecker deltas in spin 
and isospin 8 . This interaction model contains the characteristic intermediate- 
range attraction and short-range repulsion of the nucleon-nucleon force through 
Vhad and allows to include the fermionic nature of nucleons. 

The parameter set employed is displayed in tabled adjusted to reproduce 
the saturation density and binding energy per nucleon of symmetric nuclear 
matter and neutron matter, and the binding energy of finite nuclei at T = 1 
MeV. 
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Table 1: Model parameters. 
a (MeV) b (MeV) c(MeV) rf(MeV) g0 (fm) Po (MeV/c) A (fm 2 ) 



133 -47 11 29 3 120 1.5 



According to the hamiltonian eq. Q the equations of motion for each nucleon 
yield 
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We then study the solution of the equations of motion using different meth- 
ods, like the exponential fitting exponentially fitted symplectic method jl(Jj . 
but we choose the Numerov type integrator algorithm ^T] ^2] because it allows 
the best accuracy, giving order 8 or 12 for this kind of hamiltonians. 

The energy of the thermostat system is a conserved quantity 

A P 2 1 
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An effective temperature, that fluctuates around the desired initially set 
temperature T, can be defined as 
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3 Simulation Results 

The simulations presented in this work were carried out with fixed baryonic 
number density, rib, lepton fraction, Y e — ^ and a fixed number of nucleons, A, 
initially placed in a cubic box of side L = (A/n^) 1 / 3 at random. To minimize 
finite size effects we use periodic boundary conditions. Up to 1,000 particles are 
taken in this work with typical thermalization times of order 10 5 fm/c. 
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Figure 1: Energy per nucleoli versus simulation time at T = 1 MeV, rib — 
0.016/m -3 and Y e = 0.2. From top to bottom we plot E T , E NH , K, V per 
nucleoli. 



As an example of typical low density conditions we consider rib = 0.016/ to -3 , 
which is about a tenth of normal nuclear density, a temperature T = 1 MeV 
and a typical electron fraction for neutron rich matter Y e = 0.2. 

The energy of the system is conserved, according to eq. ijTUj l. as can be seen 
in the plot of energy per particle versus time in Fig. ^ The system exhibits 
characteristic oscillations in the thermostat variables due to the value of Q, that 
it is associated with the heat capacity of the system. 

In Fig. we plot the effective temperature for a system of A = 200 particles 
at rib = 0.016/m~ 3 , Y e = 0.2 and two values of Q. The upper curve (dashed 
line) corresponds to Q = 10 6 MeV (fm/c) 2 and the lower curve (solid line) to 
Q = 10 8 MeV (fm/c) 2 . The temperature is set to T=l MeV for both curves. 
The upper curve has been biased by adding an offset of 0.5 MeV for the sake 
of clarity. By decreasing Q the temperature control is better but leads to rapid 
oscillations in the energy, which must be carefully considered when studying the 
dynamical response in energy modes of the system jT3| . 

In the same way in Fig. [31 we plot the thermostat variable In s for the same 
cases as in Fig. [2] and again the dashed line corresponds to Q = 10 6 MeV 
(fm/c) 2 and the solid line to Q = 10 8 MeV (fm/c) 2 . This illustrates the slow 
time variation of Et for systems where Q is set bigger. 

In the numerical simulation the energy of the system may suffer a small drift 
with time due to the accuracy of the algorithm used. We have checked that the 
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Figure 2: Effective temperature for a configuration of 200 particles at T e ff = 
1 MeV, baryon density n b = 0.016/m" 3 and Y e = 0.2 for Q = 10 6 MeV(.fm/c) 2 
and Q = 10 8 MeV(fm/c) 2 . See text for details. 
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Figure 3: Thermostat variable Ins versus simulation time at T e ft = 1 MeV, 
baryon density rib — 0.016/m -3 and Y e — 0.2 for Q = 10 6 MeV(fm/c) 2 and 
Q = 10 s MeV(fm/c) 2 . See text for details. 
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maximum relative extended energy error at time t, denned as 
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is ^P- = 10~ 4 for time length At = 5.10 3 fm/c using timesteps dt = 0.025 fm/c 
and Q — 10 6 MeV(fm/c) 2 . This error increases somewhat with timestep size. 
We use timesteps in the range dt = 0.01 — 0.1 fm/c for our simulations in this 
work. 

4 Conclusions 

In this work we have employed molecular dynamics techniques to solve a hamil- 
tonian model with an interaction potential depending on both positions and 
momenta. We have simulated neutron rich matter at a given density and kept 
fixed the temperature by using the Nose-Hoover thermostat and solved the equa- 
tions of motion using different integration methods. We find that a Numerov 
type algorithm allows the best accuracy for this kind of hamiltonians giving 
order 8 or 12. 

We find that at the density of a tenth of nuclear saturation density Y e = 0.2 
and T=l MeV a clustered phase is formed. By changing the heat capacity of 
the system in the range Q = 10 6 — 10 8 MeV(fm/c) 2 thermalization of the 
system is achieved in times of order 10 5 fm/c. The lower the heat capacity, Q, 
the better the temperature control is. Induced oscillations in the thermostat 
variables must be considered with further detail as they could jeopardize low 
energy excitation modes. 
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